Dissipation-induced d-Wave Pairing of Fermionic Atoms in an Optical Lattice 
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We show how dissipative dynamics can give rise to pairing for two-component fermions on a 
lattice. In particular, we construct a "parent" Liouvillian operator so that a BCS-type state of 
a given symmetry, e.g. a d-wave state, is reached for arbitrary initial states in the absence of 
conservative forces. The system-bath couplings describe single-particle, number conserving and 
quasi-local processes. The pairing mechanism crucially relies on Fermi statistics. We show how 
such Liouvillians can be realized via reservoir engineering with cold atoms representing a driven 
dissipative dynamics. 



Pairing in condensed matter physics in general, and in 
atomic quantum gases in particular, is associated with 
conservative forces between particles, e.g., in Cooper 
pairs or molecular BEC pairs pQ . Lattice dynamics gives 
rise to exotic forms of pairing, such as the expected for- 
mation of d-wave Cooper pairs of fermions for a 2D Hub- 
bard model for repulsive interactions, as discussed in the 
context of high-T c superconductivity [2J, but also con- 
densates of 77-pairs [3J, and the formation of repulsively 
bound atom pairs 0]. Here we show that purely dissipa- 
tive dynamics, induced by coupling the system to a bath, 
can give rise to pairing, even in the complete absence of 
conservative forces. This "dissipative pairing" crucially 
relies on Fermi statistics and is in contrast to pairing 
arising from bath- mediated interactions (e.g., phonon- 
mediated Cooper pairing). We will discuss how reservoir 
engineering provides opportunities for experimental re- 
alisation of this dissipative pairing mechanism with cold 
atomic fermions in optical lattices [S]. 

Below we treat the example of a d-wave-paired BCS 
state of two-component fermions in two dimensions (2D), 
showing how the pairing can be generated via purely dis- 
sipative processes. A BCS-type state is the conceptually 
simplest many-body wave function describing a conden- 
sate of N paired spin-1/2 fermionic particles, |BCStv) ~ 
(dft) N / 2 \vac). On a square lattice, and assuming sin- 
glet pairs with zero center-of-mass momentum, we have 



rft 



Eq ^qC q ,T C -q4 Or S = J2i,j <Pij4,t C ll> Where C <U 

(q CT ) denotes the creation operator for fermions with 
quasimomentum q (on lattice site i) and spin o =j, |, and 
Vq (fij) the momentum (position) wave function of the 
pairs. For d-wave pairing, the pair wave function obeys 

<Pq*,q y = -<P-q v ,q* = <P-Qx,-q v ' and bel ° W We choose 

<y3 q = cosq x - cosq y or ip i;j = | Y,\= x ,y P*( 5 i,j+*x + 
6i,j- ex ) with p x = —p y = 1 corresponding to the limit of 
well localized pairs (see Fig. la), and e^ the unit lattice 
vector in A = x,y direction. For reference below we re- 
mark that in BCS theory, with pairing induced by coher- 
ent interactions, the corresponding energy gap function 
would be A q = A (cos q x — cos q y ) in the molecular limit. 
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FIG. 1: (a) Symmetry in the d-wave state, represented by 
a single offsite fermion pair exhibiting the characteristic sign 
change under spatial rotations. In a d-wave BCS state, this 
pair is delocalized over the whole lattice. (b,c) The dissipative 
pairing mechanism builds on (b) Pauli blocking and (c) der- 
ealization via phase locking, (b) Illustration of the action of 
Lindblad operators using Pauli blocking for a Neel state (see 
text), (c) The d-wave state may be seen as a derealization 
of these pairs away from half filling (shown is a cut along one 
lattice axis). 



The dissipative pairing mechanism is readily generalized 
to other pairing symmetries, such as e.g. p x +ip y [6], as 
long as the pairing is not onsite. 

While in the standard scenario BCS-type states are 
typically used as variational mean-field wavefunctions to 
describe pairing due to interactions, here the system is 
dissipatively driven towards the (pure) many-body BCS 

state, p{i) = e ct p(0) ^ |BCSjv}(BCSjv|, beginning 
from an arbitrary initial mixed state p(0). The dynam- 
ics of the density matrix for the TV-particle system p(t) 
is generated by a Liouville operator with the structure 
Cp = —\H c gp + ip-ffjfj + nj^ijepjj with non-hermitian 
effective Hamiltonian H c ^ = H — h^TliJiJi- Here, 
{j^} are non-hermitian Lindblad operators reflecting the 
system-bath coupling with strength characterized by the 
rate k. The Hamiltonian H generates unitary evolution, 
and will be set to zero for most of the discussion. The 
pure paired BCS state being the unique steady state of 
the dissipative dynamics results from the possibility to 
identify a set of operators with j^|BCSjv) = V^ [Z1H]. 
Below, we will identify these operators ji for the d- 



wave paired BCS states, and in addition study the dy- 
namics close to the final steady state, i.e., near |BCSat). 
We can then investigate the complex excitation spectrum 
of £, where, remarkably, we find a dissipative "BCS gap" 
that implies exponential approach to the steady state. 

We can readily check that the Lindblad operators jg 
generating the d-wave BCS state are given by 



j? = E pM 



+e A 



> Q c„ 



(1) 



with 2-spinor c 2 ; = {ci : \, c^) T and a a Pauli matrices with 
a = ±,z or a — x,y,z. An explicit construction is given 
below. Remarkably, these Lindblad operators, which 
generate pairing dissipatively, are bilinear and number 
conserving, thus acting on a single-particle only. They 
are also quasi-local operators, involving only a plaquette 
of nearest neighbor sites (see Fig. la). 

Before entering the more technical discussion of obtain- 
ing these J", we discuss the dynamics for states close to 
the final state |BCSat), where the physics is particularly 
transparent and analogies to the usual case of interaction- 
induced pairing in BCS theory can be made. For states 
close to |BCSjv) we can linearize the master equation dy- 
namics using a Bogoliubov-type approach. Here we take 
advantage of the fact that we know the steady state for 
our problem exactly; this knowledge can be used to con- 
struct a quadratic theory for the fluctuations on top of it. 
For this purpose it is technically convenient to give up ex- 
act particle number conservation, and to work with fixed 
phase coherent states |BCSe) = A/" -1 ^ 2 exp(e ie d^)|vac) 
instead of the number states |BCSjv) PQ, where JV = 
rici(l + Vq) ensures the normalization. The density ma- 
trix for these states, describing the dark steady state, 
factorises in momentum space since exp(e tf^)|vac) = 
]1 (l + e l6 V q c q1 x_ | )|vac). At late times, we can there- 
fore expand the state around |BCSg) by making the fac- 
torized ansatz p — FJ p q , where p q contains the modes 
±(q, a) necessary to describe pairing. We can then uti- 
lize the projection prescription p q = tr^ q p to find the 
equations of motion for the single pair density matrices 
p q in the presence of nonzero mean fields. These result 
from the coupling to other momentum modes, and their 
values are dictated by the final state properties. The 
resulting effective Hamiltonian is quadratic: 



HeB = -f E {™( C q, CT C q,<T + IVqPcq^C^) 
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with s-f = —l,si = 1 and dimensionless "gap function" 
A q = Aifq, and where the diagonal and off diagonal 



mean fields evaluate ton= | A| = 2 f (2 5 2 j 
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on the d-wave state, where the integration is over the 
Brillouin zone. We diagonalize H c s in the second line, 



introducing quasiparticle Lindblad operators 

7q,<x = (1 + <y3 q r 1/2 (C-q 1<T + S CT <y9 q Cq,-<x)- 

In this basis, the resulting master equation reads dtp = 
-iH e sp+ipHl s +J2q :C r K q 7q,(rP7q,<T- The linearized Lind- 
blad operators have analogous properties to quasiparticle 
operators familiar from interaction pairing problems: (i) 
They annihilate the (unique) steady state 7 q ,o-|BCSg) = 
0; (ii) they obey the Dirac algebra {7 q , CT ,7 q ' a >} — 
^q.q'^CT.o-' ari d zero otherwise |13j : and (iii) therefore are 
related to the original fermions via a canonical transfor- 
mation. The imaginary spectrum of the effective Hamil- 
tonian features a "dissipative pairing gap" 

ft q = K h (1 + if ) > k h. 

The dissipative gap implies an exponential approach to 
the steady d-wave BCS state for long times. This can 
be most easily seen in a quantum trajectory represen- 
tation of the master equation, where the system's time 
evolution is described by a stochastic wavefunction \ip(t)} 
evolving under a non-hermitian Hamiltonian \tp(t)) = 
e - lH ^ t \ip(0))/\\...\\ interrupted with rate k \\je\ip(t))\\ 2 
by quantum jumps \ip(t)} — > ji\ip{t))/ ||. . .|| so that p(t) = 
(|V'(t))('0(i)|)stoch (see, e.g., [H])- We thus see that (i) the 
BCS state is a "dark state" of the dissipative dynam- 
ics in the sense that j^|BCSjv) = implies that there 
will never be quantum jump, i.e. the state remains in 
|BCSjv), and (ii) states near |BCSj\r) show an exponen- 
tial decay according to the dissipative gap. Note that 
it is in marked contrast to dissipative preparation of a 
non-interacting BEC state in bosonic systems, where an 
approach polynomial in time is expected [7]. 

This convergence to a unique pure state is illustrated in 
Fig. [2] using numerical simulations for small systems. In 
Fig. 2a we show the entropy of the full density matrix for 
a small ID system as a function of time, and in Fig. 2b the 
fidelity of the BCS state for a small 2D grid, computed 
via the quantum trajectories method. 

Lindblad operators for d-wave states - We now turn 
to the construction of the Lindblad operators for the d- 
wave BCS state as given in Eq. (fTl). We will perform this 
construction first for an antiferromagnetic Neel state at 
half filling, and then generalize to the BCS state. Our 
task can be formulated as finding for a given many-body 
state |d) a set of (non-hermitian) Lindblad operators je 
so that it becomes the unique dark state, je\d) — VZ. 
Both the Neel and the BCS state have product form, 
|d) = rim^ml vac )- Thus, we note as a sufficient dark 
state condition [ji^d]^ — 0. 

There are two antiferromagnetic Neel states at 
half filling |N+) = rL e A4e.,t c !> ac >> l N -> = 



ITieA H+e m I cj ^ I vac) with A a sublattice in a two- 
dimensional bipartite (square) lattice, which differ by 
an overall spin flip. Introducing "Neel unit cell op- 
erators" S? v = c\ +ei a a c\, a = ±,e„ = {ie^ie^}, 




(a = (a, z)) is fulfilled by 



J? = YVj, 6 
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FIG. 2: Numerical illustration of the uniqueness of the steady 
state, showing evolution under the master equation with Lind- 
blad operators from Eq. |n). (a) Entropy computed exactly 
for four atoms on a 4x1 lattice, showing exponential conver- 
gence from a completely mixed state to a pure state, (b) 
Fidelity to the d-wave BCS state, (BCSivjpjBCSiv) with 4 
atoms on a 4x4 grid, computed via a quantum trajectories 
method (see text). Dashed lines show sampling error. 



whose usefulness will become apparent soon, the state 
can be written in eight different forms, |N±) = 
ILgA^lvac) = (-l) M / 2 rL 6 B^->ac), with M the 
lattice size. We then see that the Lindblad operators 
must obey [jf v , <§? J = for all i,j located on the same 
sublattice A or sublattice B, which is fulfilled for the set 



J, 



i e AotB. 



(3) 



Note that these operators can be obtained from Sf v by 
a particle-hole transformation c\ a — > ci M on the central 
site i. For the action of the operators jf v the assump- 
tion of fermionic statistics is essential, as illustrated in 
Fig. lb: they generate spin flipping transport according 



with jf = c\ +e ^a z Ci, establishing Eq. ([I]). Similar to 
above, each Jf is obtained from Df by a particle-hole 
transformation on the central site i. In fact, for these 
operators the stronger quasi-local commutation proper- 
ties with the molecular d-wave pairs holds due to Eq. 
(|]): [Jf, D|] = for all i,j,[Jf, D)] = for all i,j in the 
same sublattice, which relies again on fermionic statis- 
tics. In contrast, the operators Jf only commute with the 
symmetric superposition of all d-wave pairs D a ,. These 
operators establish coherence via phase locking between 
adjacent cloverleaves of sites. 

The dark state uniqueness for the Lindblad opera- 
tors (fTl) is equivalent to the uniqueness of the ground 
state of the associated hermitian Hamiltonian H = 
VY, l . a =±. z J f J ? for V > 0. We note that our BCS 
state shares the symmetries of H of global phase and 
spin rotations, and translation invariance. Based on the 
reasonable assumption that no other symmetries exist, 
we then expect the ground state to be unique. Note, 
however, the necessity of the full set { Jf}: Omitting e.g. 
{ Jf } gives rise to an additional discrete symmetry in H 
resulting in ground state degeneracy. These results are 
confirmed with numerical simulations for small systems 
and periodic boundary conditions, as shown in Fig. [2] . 

The above construction method allow us to find "par- 
ent" Lindblad operators for a much wider class of BCS- 
type states. For example, for a p x + ip y -wave state of 



to e.g. f. , ,. = c] , _ + Ci I , which is not possible when the . , „ . + „ + + 

o Jim j+e„,t >■<-'_ - spmless fermions, generated by p 1 ~ V, ,, p„c + „ c- with 

p x = —p- x = —ipy = ip-y = 1, the Lindblad opera- 
More generally, they can be 



antiferromagnetic order is already present. The proof of 
uniqueness of the Neel steady state up to double degen- 
eracy is then trivial: The steady state must fulfill the 
quasilocal condition that for any site occupied by a cer- 
tain spin, its neighboring sites be filled by opposite spins. 
For half filling, the only states with this property are 
|N±). The residual degeneracy can be lifted by adding a 
single operator ji = cj +e (1 + <J z )ci at arbitrary i. 

To find the Lindblad operators for the d-wave BCS 
state, we apply a similar strategy. We first rewrite the 
d-wave generator using the operators Sf 7 

i i 

v 

where p± x = —p± y — 1, and the quasilocal d-wave pair 
Df may be seen as the "d-wave unit cell operators" . Note 
the freedom of choosing a = ± in writing the state. This 
form makes the physical picture of a d-wave supcrfluid as 
delocalized antiferromagnetic order away from half filling 
[2] particularly apparent. The condition [J", ^ Dj] = 



i+e v 



tors are J t = Y^ v PvC. 
obtained for any fixed number pairing state with bilo- 
cal pairing [11]. Note, however, that the construction is 
not applicable for the onsite (singlet) pairing states - the 
analogs of Eq. |l]) become local, such that the lattice 
sites decouple and no phase coherence can be built up. 

Physical Implementation - The quasilocal and number- 
conserving form of Jf raises the possibility to realise dis- 
sipative pairing via reservoir engineering with cold atoms. 
We illustrate this, considering alkaline earth-like atoms 
[TO] with nuclear spin (e.g., / = 1/2 for m Yb), and a 
metastable 3 Po manifold which can be trapped indepen- 
dently to the ground ^o manifold. In this setting, one 
can construct a stroboscopic implementation, where the 
action of each Jf is realised successively. For clarity, we 
present this initially in ID, and choose the example of 
Jj = (cL_ i * + c]Li ^) c i- The implementation is depicted 
in Fig. 3: (i) The 3 Po state is trapped in a lattice of three 
times the period as that for the 1 So state, defining blocks 
of three sites in the 1 So lattice. Using this, any J, atom in 
1 So on the central site is excited to the f state of the 3 Po 
manifold, (ii) By adding an additional potential the traps 




FIG. 3: (a) Level scheme for alkaline earth atoms with I = 
1/2, showing excitation to a metastable level whilst flipping 
the nuclear spin, and induced decay by coupling to the x Pi 
level, b) Illustration of J+ implementation in ID: (i) A longer 
period lattice for 3 Po identifies a triple of wells, and atoms 
from the central level are transferred to the 3 Po manifold 
with spin flip, (ii) The 3 Po potential wells are adiabatically 
split into two; (iii) Decay is induced, returning the atom to 
the So level via coupling to a lossy cavity mode. 



For k — > and interaction V > this Hamiltonian can 
be identified as a parent Hamiltonian [11] with |BCSjv) 
as unique stable ground state and gapped positive defi- 
nite excitation spectrum. This parent Hamiltonian could 
be realised via a similar procedure to the induced dissi- 
pation, replacing the decay in step (iii) by induced in- 
teractions between atoms. This opens the possibility to 
use the d-wave state as an initial state for the prepa- 
ration of the ground state of the Fermi-Hubbard model 
by a suitable adiabatic passage [12]. Here, one can take 
advantage of the fact that (i) in the initial stages the sys- 
tem is protected by a gap ~ 0.72]/, and (ii) the d-wave 
state has identical symmetry and similar energy to the 
conjectured Fermi-Hubbard ground state away from half 
filling. Thus, since no phase transition has to be crossed, 
a d-wave superfhiid gap protection persists through the 
whole passage path. 

We thank A. Gorshkov, K. Hammerer, B. Kraus and 
A. Kantian for discussions. This work was supported by 
the Austrian Science Fund through SFB F40 FOQUS and 
EUROQUAM_DQS (I118-N16), and the EU through IP 
AQUTE. 



for 3 Po are divided so that atoms confined in them over- 
lap the right and left sites of the 3-site block for 1 S . (iii) 
Dissipation is induced via spontaneous decay, obtained 
by coupling atoms in the 3 Po state off-resonantly to the 
1 Pi state, as depicted in Fig. 3a, with coupling strength 
O, and detuning A. If we couple the 1 Srj- 1 Pi transition 
to a cavity mode with linewidth T and vacuum Rabi fre- 
quency g, then the decay will be coherent over the triple 
of sites. In the limit A 3> il and r ~^> -&■ , we obtain an 

effective decay rate T G ff = -^S^ ~ 9kHz for typical pa- 
rameters, which bounds the effective dissipative rate for 
the stroboscopic process, k. Provided atoms remain in 
the lowest band, Fermi statistics will be respected, and 
coherent dynamics during this process can be neglected 
in a deep lattice for small scattering lengths. 

This operation can occur in parallel for different 3- 
site blocks, and should be repeated with the superlattice 
shifted for other central sites. Similar operations com- 
bined with rotations of the nuclear spin before and after 
these operations allows implementation of J~ and Jf. In 
2D 3x3 plaquettes are defined by the appropriate super- 
lattice potential for the 3 Po level, and the adiabatic ma- 
nipulation of the potential in step (ii), should be adjusted 
to ensure that the correct relative phases are obtained for 
atoms transported in orthogonal directions. 

The d-wave parent Hamiltonian — As a final remark, 
we note that the effective Hamiltonian above can be gen- 
eralized to include a coherent interaction V, 



H c , 



(v- 



^)E J - 
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